#include <assert.h>
#include <stdio.h>

#include "../src/mstep.c"

static void
fill (double*           x,
      double          (*func) (void),
      unsigned int      n)
{
  for (unsigned int i = 0; i < n; ++i)
    { 
      x[i] = func ();
    }
}

static void
fill_r (Rating*         x,
        unsigned int    n_workers,
        unsigned int    n_labels,
        unsigned int    n_images,
        unsigned int    n)
{
  for (unsigned int i = 0; i < n; ++i)
    { 
      x[i].worker = (lrand48 () >> 16) % n_workers;
      x[i].image = (lrand48 () >> 16) % n_images;
      x[i].rating = (lrand48 () >> 16) % n_labels;
    }
}

static void
test_to_from_identity_once (void)
{ 
  /* (kappa, logalpha, tau, gamma) = 
     from_x (kappa, logalpha, tau, gamma, to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];

  double newkappa;
  double newlogalpha[n_images];
  double newtau[n_workers * (n_labels - 1)];
  double newgamma[n_labels - 1];
  double newrho;
  double newloglambda[n_workers];

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);

  newkappa = drand48 ();
  fill (newlogalpha, drand48, n_images);
  fill (newtau, drand48, n_workers * (n_labels - 1));
  fill (newgamma, drand48, n_labels - 1);
  newrho = drand48 ();
  fill (newloglambda, drand48, n_workers);

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);
  from_x (&newkappa, newlogalpha, n_workers, n_labels, newtau, n_images, newgamma, &newrho, newloglambda, x);
          
  gsl_vector_free (x);

  assert (kappa == newkappa);
  assert (memcmp (logalpha, newlogalpha, n_images * sizeof (logalpha[0])) == 0);
  assert (memcmp (tau,
                  newtau,
                  n_workers * (n_labels - 1) * sizeof (tau[0])) == 0);
  assert (memcmp (gamma, newgamma, (n_labels - 1) * sizeof (gamma[0])) == 0);
  assert (rho == newrho);
  assert (memcmp (loglambda, newloglambda, n_workers * sizeof (loglambda[0])) == 0);
}

static void
test_to_from_inverse (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_to_from_identity_once ();
    }
}

static void
test_qfunc_f_once (void)
{
  /* qfunc (kappa, logalpha, tau, gamma) = 
     -qfunc_f (to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  double hq = 0;

  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .priorkappa = &priorkappa.base,
                               .priorlogalpha = &priorlogalpha.base,
                               .priortau = &priortau.base,
                               .priorgamma = priorgamma,
                               .priorrho = &priorrho.base,
                               .priorloglambda = &priorloglambda.base,
                               .hq = hq };

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);

  double q = qfunc (kappa, logalpha, n_workers, n_labels,
                    tau, n_images, gamma, rho, loglambda, pz, 
                    ratings, n_ratings, &priorkappa.base,
                    &priorlogalpha.base, &priortau.base,
                    priorgamma, &priorrho.base, &priorloglambda.base, hq);

  double qf = qfunc_f (x, &qfunc_params);

  assert (fabs (q + qf) <= 1e-8 * (1 + fabs (q) + fabs (qf)) ||
          (fprintf (stderr, "%g ?= %g\n", q, qf), 0));

  gsl_vector_free (x);
}

static void
test_qfunc_f (void)
{ 
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_qfunc_f_once ();
    }
}

static void
test_qfunc_df_once (void)
{
  /* dqfunc (kappa, logalpha, tau, gamma) = 
     -qfunc_df (to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double dkappa;
  double dlogalpha[n_images];
  double dtau[n_workers * (n_labels - 1)];
  double dgamma[n_labels - 1];
  double drho;
  double dloglambda[n_workers];
  double newdkappa;
  double newdlogalpha[n_images];
  double newdtau[n_workers * (n_labels - 1)];
  double newdgamma[n_labels - 1];
  double newdrho;
  double newdloglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());
    
  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  dkappa = drand48 ();
  fill (dlogalpha, drand48, n_images);
  fill (dtau, drand48, n_workers * (n_labels - 1));
  fill (dgamma, drand48, n_labels - 1);
  drho = drand48 ();
  fill (dloglambda, drand48, n_workers);
  newdkappa = drand48 ();
  fill (newdlogalpha, drand48, n_images);
  fill (newdtau, drand48, n_workers * (n_labels - 1));
  fill (newdgamma, drand48, n_labels - 1);
  newdrho = drand48 ();
  fill (newdloglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .priorkappa = &priorkappa.base,
                               .priorlogalpha = &priorlogalpha.base,
                               .priortau = &priortau.base,
                               .priorgamma = priorgamma,
                               .priorrho = &priorrho.base,
                               .priorloglambda = &priorloglambda.base };

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);
  gsl_vector* g = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);

  dqfunc (&dkappa, dlogalpha, dtau, dgamma, &drho, dloglambda,
          kappa, logalpha, n_workers, n_labels,
          tau, n_images, gamma, rho, loglambda, 
          pz, ratings, n_ratings, &priorkappa.base,
          &priorlogalpha.base, &priortau.base,
          priorgamma, &priorrho.base, &priorloglambda.base);

  qfunc_df (x, &qfunc_params, g);

  from_x (&newdkappa, newdlogalpha, n_workers, n_labels, newdtau, n_images, newdgamma, &newdrho, newdloglambda, g);

  assert (fabs (dkappa + newdkappa) <= 
          1e-8 * (1 + fabs (dkappa) + fabs (newdkappa)) ||
          (fprintf (stderr, "%g ?= %g\n", dkappa, newdkappa), 0));

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (fabs (dlogalpha[i] + newdlogalpha[i]) <=
              1e-8 * (1 + fabs (dlogalpha[i]) + fabs (newdlogalpha[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dlogalpha[i], newdlogalpha[i]), 0));
    }

  for (unsigned int i = 0; i < n_workers * (n_labels - 1); ++i)
    {
      assert (fabs (dtau[i] + newdtau[i]) <=
              1e-8 * (1 + fabs (dtau[i]) + fabs (newdtau[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dtau[i], newdtau[i]), 0));
    }

  for (unsigned int i = 0; i < (n_labels - 1); ++i)
    {
      assert (fabs (dgamma[i] + newdgamma[i]) <=
              1e-8 * (1 + fabs (dgamma[i]) + fabs (newdgamma[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dgamma[i], newdgamma[i]), 0));
    }

  assert (fabs (drho + newdrho) <= 
          1e-8 * (1 + fabs (drho) + fabs (newdrho)) ||
          (fprintf (stderr, "%g ?= %g\n", drho, newdrho), 0));

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      assert (fabs (dloglambda[i] + newdloglambda[i]) <=
              1e-8 * (1 + fabs (dloglambda[i]) + fabs (newdloglambda[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dloglambda[i], newdloglambda[i]), 0));
    }

  gsl_vector_free (g);
  gsl_vector_free (x);
}

static void
test_qfunc_df (void)
{ 
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_qfunc_df_once ();
    }
}

static void
test_mstep_no_data_once (void)
{
  /* without data mstep should converge to prior maximum */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  double hq = 0;

  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  mstep (&kappa, logalpha, n_workers, n_labels,
         tau, n_images, gamma, &rho, loglambda,
         pz, priorz, NULL, 0, 
         &priorkappa.base, &priorlogalpha.base,
         &priortau.base, priorgamma, &priorrho.base, &priorloglambda.base,
         hq, 100, 1e-8);

  assert (   fabs (kappa - priorkappa.mean) 
          <= 1e-3 * (1 + fabs (kappa) + fabs (priorkappa.mean))
          || (fprintf (stderr, "%g ?= %g\n", kappa, priorkappa.mean), 0));

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logalpha[i] - priorkappa.mean)
              <= 1e-3 * (1 + fabs (logalpha[i]) + fabs (priorkappa.mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, logalpha[i], priorkappa.mean), 0));
    }

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      assert (   fabs (gamma[i] - normalpriorgamma[i].mean)
              <= 1e-3 * (1 + fabs (gamma[i]) + fabs (normalpriorgamma[i].mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], normalpriorgamma[i].mean), 0));
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j + 1 < n_labels; ++j)
        {
          assert (   fabs (tau[i * (n_labels - 1) + j] - normalpriorgamma[j].mean)
                  <= 1e-3 * (1 + fabs (tau[i * (n_labels - 1) + j]) + fabs (normalpriorgamma[j].mean))
                  || (fprintf (stderr, "%u: %g ?= %g\n", i, tau[i * (n_labels - 1) + j], normalpriorgamma[j].mean), 0));
        }
    }

  assert (   fabs (rho - priorrho.mean) 
          <= 1e-3 * (1 + fabs (rho) + fabs (priorrho.mean))
          || (fprintf (stderr, "%g ?= %g\n", rho, priorrho.mean), 0));

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      assert (   fabs (loglambda[i] - priorrho.mean)
              <= 1e-3 * (1 + fabs (loglambda[i]) + fabs (priorrho.mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, loglambda[i], priorrho.mean), 0));
    }
}

static void
test_mstep_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_no_data_once ();
    }
}

typedef struct _MstepTestDatum MstepTestDatum;
struct _MstepTestDatum 
{
  struct
    {
      double kappa;
      double logalpha[3];
      double tau[4*(5-1)];
      double gamma[(5-1)];
      double rho;
      double loglambda[4];
      double pz[3*5];
      double priorz[5];
      Rating ratings[8];
      double priorlogalphastddev;
      double priortaustddev;
      double priorgamma1mean;
      double priorgamma1stddev;
      double priorgamma2mean;
      double priorgamma2stddev;
      double priorgamma3mean;
      double priorgamma3stddev;
      double priorgamma4mean;
      double priorgamma4stddev;
      double priorkappamean;
      double priorkappastddev;
      double priorrhomean;
      double priorrhostddev;
      double priorloglambdastddev;
    }                           input;

  struct
    {
      double q;
      struct
        {
          double kappa;
          double logalpha[3];
          double tau[4*(5-1)];
          double gamma[(5-1)];
          double rho;
          double loglambda[4];
        }    v;
    }                           output;
};

static void
test_mstep (void)
{
  MstepTestDatum data[] = 
    {{{0.353235, {0.580682, 0.239041, 0.703317}, {0.874543, 0.39361,0.26239, 0.692018, 0.0578394, 0.823095, 0.652613, 0.350876,0.378499, 0.71082, 0.00135051, 0.564658, 0.822995, 0.0477758,0.915056, 0.0492184}, {0.174056, 0.858161, 0.407386, 0.23479},0.820821, {0.277479, 0.168345, 0.531473, 0.946278}, {0.883869,0.905955, 0.839455, 0.888439, 0.060774, 0.253342, 0.488579,0.50994, 0.349954, 0.251992, 0.923921, 0.686945, 0.302178,0.336936, 0.874703}, {0.512889, 0.444017, 0.92955, 0.639913,0.692068}, {{1, 2, 0}, {2, 0, 3}, {2, 2, 0}, {0, 1, 2}, {0, 0,2}, {0, 0, 4}, {3, 2, 0}, {1, 1, 1}}, 0.416538, 1.0112, 0.10844,0.99579, 0.282669, 1.10525, 0.268985, 1.10735, 0.221896, 0.851908,0.780406, 0.597412, 0.871942, 0.599916,1.10648}, {-66.547968318629761662, {-0.38828329607898211253,{-0.33704417964502518813, -0.53076901907036623860,-0.86518449774779550928}, {0.10569768636912235433,-0.59023070716028070204, 1.6580108914510073784,0.29034446572714102268, 0.59101419816813500073,0.40975087331111164761, 1.3482744419030647858,0.54800288844031154633, 0.58614330682716366001,0.12156855819754325191, 0.96835182403226451612,0.67659671384334784869, 0.67086325299201850684,0.30592450134432813016, 1.2765762711236082118,0.50987580255579623557}, {0.41054649853412294553,0.099983214991049796510, 1.1327357047284318026,0.43214691200740774383}, -0.30656475872292192382,{-0.66395436106059783573, -1.6586007403601577249,-1.4476664989302337672, -1.4650963723544891626}}}}, {{0.660467,{0.569763, 0.0129796, 0.981782}, {0.147578, 0.125746, 0.0834297,0.34187, 0.45551, 0.959208, 0.322225, 0.23343, 0.70972, 0.676538,0.466975, 0.964445, 0.852368, 0.454643, 0.865067,0.184039}, {0.504956, 0.582701, 0.515152, 0.327554},0.844489, {0.0129377, 0.502172, 0.345772, 0.696911}, {0.887192,0.418742, 0.00390203, 0.241401, 0.927984, 0.0965175, 0.770472,0.531681, 0.251446, 0.629543, 0.806027, 0.679313, 0.796803,0.764475, 0.621988}, {0.174357, 0.214102, 0.249324, 0.294434,0.329868}, {{2, 1, 2}, {1, 0, 0}, {2, 2, 2}, {2, 2, 3}, {0, 1,2}, {2, 1, 1}, {1, 0, 2}, {2, 1, 2}}, 0.451164, 0.997152,0.948663, 0.882958, 0.313973, 0.578409, 0.94476, 0.641557,0.385989, 0.481892, 0.174289, 1.10988, 0.134543, 0.852349,0.618262}, {-53.876800700824913364, {-1.3153203387660448968,{-1.4594678441123200864, -1.4447704108783462642,-1.2878691410576526638}, {0.82006202411067869643,0.15073870769944347986, 1.3135797448649448854,0.57925245624673608052, 0.99288561281722748354,0.29002480130515693046, 1.3818857137173358389,0.61407480863573753792,0.63389775086020222162, -0.026076314077639740337,1.4827525496557898441, 0.88848267545448498887,0.85530442451089070763, 0.22568917096287833093,1.1929137226052061053,0.51285950976716705334}, {0.85530442445374730310,0.22568917114649646708, 1.1929137226472585864,0.51285950990080837974}, -0.74399254696977262695,{-0.78176429100468986902, -1.0146887369852856503,-0.89776561393871090833, -0.74399254701727034588}}}}, {{0.180562,{0.33774, 0.837874, 0.746274}, {0.00620448, 0.123638, 0.58855,0.451839, 0.676336, 0.922473, 0.841399, 0.503177, 0.0433783,0.608501, 0.512989, 0.558416, 0.651822, 0.222512, 0.281097,0.384128}, {0.791946, 0.0879692, 0.678748, 0.0158663},0.611385, {0.750229, 0.840874, 0.269593, 0.60518}, {0.626592,0.252323, 0.817753, 0.928844, 0.704118, 0.410924, 0.314576,0.885466, 0.0956177, 0.897935, 0.75616, 0.233644, 0.873106,0.616838, 0.372032}, {0.441698, 0.785137, 0.93809, 0.356166,0.830313}, {{3, 1, 1}, {1, 1, 1}, {3, 1, 4}, {0, 2, 4}, {0, 0,1}, {2, 2, 2}, {2, 1, 3}, {2, 2, 1}}, 0.284907, 0.347217,0.0865733, 0.475133, 0.408315, 1.09489, 0.26882, 0.546289,0.704197, 0.683969, 0.954243, 0.660823, 0.608579, 0.786034,0.448083}, {-39.321626734275088067, {-0.35548749425279570425,{-0.45788418504825998947, -0.44936272535739412281,-0.40267015223566480018}, {-0.031665333250674341893,1.3730668229903959839, 0.45318960905929071689,0.75493156470128701075, -0.021874797574626214637,1.3786641656187735473, 0.47195431396547129448,0.79087487799089250613, -0.036694792020035855890,1.3512122727177587405, 0.47876814985120584902,0.82035430014529402010, -0.028517086158621219036,1.3689915956033972524, 0.45232478225956492486,0.75719295681062201977}, {-0.015994316848792126221,1.3444476199032751840, 0.44614993208012059893,0.77619951241612859201}, -1.2444993564734760781,{-1.3716353892843347200, -1.4341692920021593118,-1.4574482966629499236, -1.3169282275437748426}}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution priorlogalpha = 
        parameterized_normal_distribution (data[i].input.priorlogalphastddev);

      ParameterizedNormalDistribution priortau =
        parameterized_normal_distribution (data[i].input.priortaustddev);

      NormalDistribution priorgammanormal[4] = {
          normal_distribution (data[i].input.priorgamma1mean,
                               data[i].input.priorgamma1stddev),
          normal_distribution (data[i].input.priorgamma2mean,
                               data[i].input.priorgamma2stddev),
          normal_distribution (data[i].input.priorgamma3mean,
                               data[i].input.priorgamma3stddev),
          normal_distribution (data[i].input.priorgamma4mean,
                               data[i].input.priorgamma4stddev)
      };

      Distribution* priorgamma[4] = {
          &priorgammanormal[0].base,
          &priorgammanormal[1].base,
          &priorgammanormal[2].base,
          &priorgammanormal[3].base
      };

      NormalDistribution priorkappa = 
        normal_distribution (data[i].input.priorkappamean,
                             data[i].input.priorkappastddev);

      NormalDistribution priorrho = 
        normal_distribution (data[i].input.priorrhomean,
                             data[i].input.priorrhostddev);

      ParameterizedNormalDistribution priorloglambda =
        parameterized_normal_distribution (data[i].input.priorloglambdastddev);

      double kappa = data[i].input.kappa;
      double logalpha[3];
      double tau[4*(5-1)];
      double gamma[(5-1)];
      double rho = data[i].input.rho;
      double loglambda[4];

      memcpy (logalpha, data[i].input.logalpha, 3 * sizeof (double));
      memcpy (tau, data[i].input.tau, 4*(5-1) * sizeof (double));
      memcpy (gamma, data[i].input.gamma, (5-1) * sizeof (double));
      memcpy (loglambda, data[i].input.loglambda, 4 * sizeof (double));

      double hq = 0;

      for (unsigned int k = 0; k < 3 * 5; ++k)
        {
          hq -= data[i].input.pz[k] * log (data[i].input.pz[k]);
        }

      double q = 0;

      for (unsigned int n = 0; n < 10; ++n)
        {
          q = mstep (&kappa,
                     logalpha,
                     4,
                     5,
                     tau,
                     3,
                     gamma,
                     &rho,
                     loglambda,
                     data[i].input.pz,
                     data[i].input.priorz,
                     data[i].input.ratings,
                     8,
                     &priorkappa.base,
                     &priorlogalpha.base,
                     &priortau.base,
                     priorgamma,
                     &priorrho.base,
                     &priorloglambda.base,
                     hq,
                     10000,
                     1e-16);
        }

      assert (      fabs (q - data[i].output.q) 
                 <= 1e-3 * (1 + fabs (q) + fabs (data[i].output.q))
              || (fprintf (stderr, "fail: %g != %g\n", q, data[i].output.q),
                  0));

      assert (   fabs (kappa - data[i].output.v.kappa)
              <= 1e-2 * (1 + fabs (kappa) + fabs (data[i].output.v.kappa))
              || (fprintf (stderr, "%g ?= %g\n", kappa, data[i].output.v.kappa), 0));

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (   fabs (logalpha[j] - data[i].output.v.logalpha[j])
                  <= 1e-2 * (1 + fabs (logalpha[j]) + fabs (data[i].output.v.logalpha[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, logalpha[j], data[i].output.v.logalpha[j]), 0));
        }

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k + 1 < 5; ++k)
            {
              assert (   fabs (tau[j * (5 - 1) + k] - data[i].output.v.tau[j * (5 - 1) + k])
                      <= 2e-2 * (1 + fabs (tau[j * (5 - 1) + k]) + fabs (data[i].output.v.tau[j * (5 - 1) + k]))
                      || (fprintf (stderr, "%u: %g ?= %g\n", j, tau[j * (5 - 1) + k], data[i].output.v.tau[j * (5 - 1) + k]), 0));
            }
        }

      for (unsigned int j = 0; j + 1 < 5; ++j)
        {
          assert (   fabs (gamma[j] - data[i].output.v.gamma[j])
                  <= 2e-2 * (1 + fabs (gamma[j]) + fabs (data[i].output.v.gamma[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, gamma[j], data[i].output.v.gamma[j]), 0));
        }

      assert (   fabs (rho - data[i].output.v.rho)
              <= 1e-2 * (1 + fabs (rho) + fabs (data[i].output.v.rho))
              || (fprintf (stderr, "%g ?= %g\n", rho, data[i].output.v.rho), 0));

      for (unsigned int j = 0; j < 4; ++j)
        {
          assert (   fabs (loglambda[j] - data[i].output.v.loglambda[j])
                  <= 1e-2 * (1 + fabs (loglambda[j]) + fabs (data[i].output.v.loglambda[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, loglambda[j], data[i].output.v.loglambda[j]), 0));
        }
    }
}

int
main (void)
{
  srand48 (69);

  test_to_from_inverse ();
  test_qfunc_f ();
  test_qfunc_df ();
  test_mstep_no_data ();
  test_mstep ();

  return 0;
}
